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Abstract 


The  nonlinear,  large  amplitude  free  vibration  of  composite  helicopter  blades  under 
large  static  deflection  is  investigated  analytically.  A  new  model  capable  of  handling 
large  amplitudes  as  well  as  large  deflections  was  developed,  based  on  the  work  in  a 
previous  report  by  Minguet.  The  model  can  deal  with  large  displacements  and  rota¬ 
tions  by  use  of  Euler  angles  and  can  account  for  structural  couplings  such  as  bending- 
twist  and  extension-twist.  The  reduction  of  this  large  deflection  model  to  a  commonly 
used  moderate  deflection  model  is  also  shown.  A  Newton-Raphson  type  iterative  so¬ 
lution  technique  based  on  numerical  integration  of  the  basic  large  deflection  equations 
is  seen  effective  for  the  present  analysis.  Two  different  lay-ups  [0/90]3,,  [45/0],,  of 
graphite/epoxy  flat  beams  have  been  selected  to  demonstrate  the  large  amplitude 
analysis.  The  behavior  of  the  first  and  second  bending,  the  first  fore-and  aft,  and  the 
first  torsional  modes  is  presented  as  tip  static  deflection  and  tip  amplitudes  increase. 
It  is  found  that  both  large  static  deflection  and  large  amplitudes  can  affect  the  fore- 
and-aft  and  torsion  modes  significantly,  but  bending  modes  sure  not  influenced  much 
by  the  geometrical  nonlinearities. 
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Chapter  1 
Introduction 


The  behavior  of  nonlinear,  large  amplitude  free  vibration  of  composite  helicopter  rotor 
blades  under  large  static  deflections  is  investigated.  A  previous  report  by  Minguet, 
Ref.  1  indicates  that  under  large  static  deflections,  natural  frequencies  and  mode 
shapes  of  the  blade,  particularly  those  of  the  fore-and-aft  and  torsion  modes  show 
interesting  trends  that  are  not  apparent  from  the  characteristics  of  undeformed 
cantilever  beams.  The  influence  of  the  large  static  deflections  on  the  modes  was 
found  by  linearizing  the  governing  equations  of  motion  around  a  given  static  position 
to  yield  the  small  amplitude  vibrations  of  the  beam  around  that  large  static  position. 

In  the  present  analysis,  the  amplitudes  of  motion  are  also  allowed  to  be  large,  and 
emphasis  is  given  on  how  the  vibrational  behavior  of  the  blades  is  affected  by  not 
only  the  static  deflection  but  also  the  amplitude  level  at  the  tip.  This  type  of  free 
vibration  analysis  should  give  insight  into  more  general  aeroelastic  problems  where 
large  amplitude  motion  is  accompanied  by  nonlinear  exciting  forces  such  as  nonlinear 
periodic  forces  due  to  gravity  or  aerodynamic  loads  due  to  dynamic  stall.  A  simple 
such  analysis,  dealing  only  with  geometrical  nonlinearities  of  the  rigid  blade,  was  given 
by  Chopra  and  Dugundji  in  Ref.  2.  More  recently  Dunn  and  Dugundji  have  given 
another  such  analysis,  this  time  dealing  only  with  aerodynamic  stall.  A  fully  nonlinear 
aeroelastic  analysis  involving  both  structural  and  aerodynamic  nonlinearities  would  be 
of  interest.  Flexible  helicopter  blades  are  good  examples  in  which  these  nonlinearities 
play  important  roles,  and  the  present  analysis  should  serve  as  an  introduction  to  the 
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understanding  of  such  complex  phenomenon. 

There  exist  two  types  of  nonlinear  helicopter  blade  equations  that  can  be  readily 
available  for  the  purpose  of  present  analysis;  the  equations  that  are  based  on  various 
geometrical  ordering  schemes,  and  the  ones  that  we  not  based  on  ordering  schemes. 
The  former  group  of  equations  approximate  large  displacements  and  rotations  mostly 
up  to  second  order(e.g.  Ref.  4,  5)  while  the  latter  group  preserve  the  complete 
nonlinearities  in  them  (e.g.  Ref.  1).  Since  strong  couplings  between  various  static  and 
dynamic  parts  of  the  equations  are  expected  in  the  nonlinear  large  amplitude 
vibrations,  the  set  of  complete  nonlinear  equations  of  the  latter  group  is  preferred. 
The  nonlinear  equations  derived  by  Minguet  in  Ref.  1  are  used  here  for  their  simplicity 
and  immediate  availability  for  analysis  of  composite  blades.  However,  to  illustrate 
correspondency  between  these  two  different  types  of  equations,  an  attempt  is  made  to 
reduce  the  nonlinear  equations  by  Minguet  t  ne  second  order  equations  for  moderate 
deflections  that  are  given  by  Hodges  and  Dowell  in  Ref.  4,  and  Boyd  in  Ref.  5. 

A  new  technique  based  on  harmonic  balance  and  iterative  Newton- Raphson 
algorithm  is  introduced  to  solve  for  the  modes  and  their  frequencies  as  functions  of 
amplitudes  of  interest  under  moderate  to  large  static  deflections.  Results  of  numerical 
analysis  are  given  for  two  lay-ups  [0/90j3,,  and  [45/0],  of  graphite/epoxy  composite 
beams  under  various  static  deflections. 

All  assumptions  made  earlier  in  Ref.  1  are  retained  throughout  the  analysis.  They 
are,  the  blade  itself  is  long  enough  to  be  treated  as  a  one-dimensional  model,  shear 
deformation  can  be  neglected,  and  warping  of  the  cross  section  of  the  blade  can 
be  neglected.  Also,  material  nonlinearity  is  not  considered  here.  As  indicated  by 
Friedmann  in  Ref.  6,  this  model  has  some  limitations  since  it  does  not  include  shear 
deformation  and  warping  of  the  blades  which  may  be  present  to  a  small  extent  in 
realistic  helicopter  blades.  However,  it  is  only  a  matter  of  refining  to  include  such 
structural  effects,  and  for  the  purposes  of  current  analysis  the  model  is  found  to 
be  adequate  to  show  the  basic  characteristics  of  large  amplitude  free  vibration  of 


composite  helicopter  blades. 


Chapter  2 

Analytic  Modeling 


2.1  Basic  Equations 

There  are  twelve  first-order,  nonlinear  differential  equations  that  describe  the  statics 
and  dynamics  of  composite  blades  completely.  For  thorough  derivation  of  the  equa¬ 
tions  see  Ref.  1.  All  the  equations  are  derived  based  on  the  following  transformation 
matrix  that  transforms  the  global  coordinate  x,y,  z  into  the  local  one  £ ,r ?,  (  (see  figure 
1),  i.  e.  , 


!h 


COS  P  cos  ip 

cos  P  sin  ip 

sin  p 

—  cos  8  sin  ip 

cos  8  cos  ip 

sin  8  cos  P 

—  sin  8  sin  P  cos  ip 

—  sin  8  sin  P  sin  ip 

sin  8  sin  ip 

—  sin  8  cos  ip 

cos  8  cos  P 

—  cos  8  sin  P  cos  ip 

—  cos  8  sin  P  sin  ip 

(2.1) 


Here  ip,P,9  are  the  local  Euler  angles.  The  transformation  matrix  is  orthogonal  and 
related  to  the  rotation  (or  curvature)  matrix  as  follows. 
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pr1  =  pir 

®-wm 


where 


3*  .  „ 

"«  =  57  +  s'“'3 


0 

Kc 

Kn 

-*c 

0 

*< 

-*C 

0 

dip 

ds 

(twist  rate) 


/j  ,  Q  o 

k-  =  —  cos  a  — — *-  sin  6  cos  p  — 


■  a  90  .  -dip 

Hr  =  sin  6  — — t-  cos  0  cos  p  — — 
03  03 


(bending  about  tj  axis) 

(2-4) 

(bending  about  £  axis) 


Inverting  the  above  differential  equation  yields 

86 

—  =  —  sin  6  tan  (3  —  cos  6  tan  (3  kq 

03 

80  ... 

—  =  —  cos  6  Kn  +  sin  6  Kr 

03 

dtp  sin  6  cos  6 

=  - n  Kn  + - a  K( 

Os  cos  p  cos  p 

The  global  displacements  are  related  to  Euler  angles  via 


— -  =  (1  4-  e)  cos 0  cosip 
Os 

=  (1  +  e)  cos/?  sin  ip  (2.6) 

dz  .  . 

—  =  (1  +  e)  sm 

where  e  is  the  axial  strain  along  the  reference  line.  In  addition  to  the  above  six 
compatibility  equations,  one  has  to  consider  equilibrium  of  forces  and  moments  of  the 
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beam.  The  equilibrium  equations  can  be  written  either  in  global  or  local  coordinates. 
Here  they  are  written  in  local  coordinate  in  order  to  take  into  account  the  large 
deformation  of  the  beam  in  space.  The  first  three  differential  equations  that  describe 
the  equilibrium  of  the  local  force  resultants  Fi,F2,  F3  are 

dF, 

-r - *£  F2  +  F3  +  Til  Pi  +  T\2  Pv  +  T13  pz  +  p\  =  0 

OS 

QF 

-3—  +  F\  —  F3  +  T21  Pi  +  T22  Pv  +  T2 3  pz  +  P2  =  0  (2.7) 

03 

dF-x 

- rtf?  F\  +  F2  +  T31  Pi  +  T32  Pv  +  T33  Pi  +  P3  =  0 

with 


p^,  :  applied  load  vector  in  local  axis  =  pi,  ?2,  P3 

pc  :  applied  load  vector  in  global  axis  =  pr,  pv,  pz 

The  other  three  differential  equations  describe  the  equilibrium  of  the  local  moment 
resultants  Mi,  M2,  M3. 

d 

M2  +  Kn  M3  +  Tnmz  +  +  Ti3m,  +  rri!  =  0 

03 

Q  1/ 

- - - h  M\  —  AC(  M3  +  T2\  HI x  +  T22  TTly  T23  TTlx  +  7712  F3  =  0  (2.8) 

03 

dMi 

— r -  —  K,)  Ml  +  M2  4-  T31  mz  +  T32  TTly  +  T33  mz  +  7713  +  F2  =  0 

03 

with 


mi  :  applied  moment  vector  in  local  axis  =  mi,  m2,  m3 
me  :  applied  moment  vector  in  global  axis  =  mz,  mv,  mt 

In  helicopter  problems,  generally  two  kinds  of  loadings  are  arised;  inertial  loads  that 
include  normal  and  angular  acceleration,  Coriolis  acceleration,  centrifugal  and  gravi¬ 
tational  forces,  and  aerodynamic  loads  that  include  both  steady  and  unsteady  parts. 
The  former  group  usually  appears  as  the  global  pc,  tog  while  the  latter  group  appears 
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as  the  local  pi,  and  mi.  In  the  present  analysis  only  the  normal,  angular  acceleration 
and  gravitational  loads  are  considered.  Hence  for  a  blade  without  mass  centroid  offset 


px  =  -mi 


Pv  =  -my 

Pj  =  -m  z  -  mg  (2.9) 

mx  =  rriy  =  m,  —  0 


and 


Pi  =  P2  —  p3  =  0 


rrit 


(2.10) 


m.2  =  m3  =  0 


Finally,  a  set  of  generalized  stress-strain  relations  are  incorperated  via  the  following 
six  linear  equations. 


r  Fi ' 

Eu  Eu  E i3  E is  E i6 

e 

f2 

Eu  Eu  Eh  Eh  Eie 

T'f’j 

F 3 

E33  E34  £35  £36 

7f( 

M\ 

£44  £45  E  46 

*£ 

M2 

SYM  £55  £5 e 

.  m3  . 

Ess 

.  . 

(2.11) 


Here  7^,7^  represent  the  two  transverse  shear  strains.  In  its  most  general  case,  the 
above  stiffness  matrix  can  be  full,  i.  e.  there  can  be  couplings  between  all  of  three 
force  resultants,  three  moment  resultants  and  all  of  six  strain  components.  However  in 
consistancy  with  the  earlier  assumptions  of  a  Bernoulli- Euler  beam,  the  calculations 
of  the  two  shear  strains  are  completely  ignored  during  the  current  analysis. 


2.2  Reduction  of  Basic  Equations  for  Moderate 
Deflections 

Before  proceeding  with  the  large  amplitude  vibration  solution  of  the  basic  nonlin¬ 
ear  equations  presented  in  the  previous  section,  the  equations  of  motion  in  u,  v ,  w, 
and  <f>  that  were  derived  by  Hodges  and  Dowell,  Ref.  4,  and  Boyd,  Ref.  5,  for  mod¬ 
erate  deflections  will  be  rederived  from  the  twelve  general  nonlinear  equations  2.5 
through  2.8.  Only  the  case  of  isotropic  blade  with  no  mass  centroid  offset  is  consid¬ 
ered  here  for  illustration.  In  this  way,  the  approximations  of  the  moderate  deflection 
analysis  can  be  assessed. 

The  first  step  in  the  reduction  process  is  to  rewrite  the  force  and  moment  equi¬ 
librium  equations  in  global  x,  y,  z  directions  instead  of  local  £,  q,  £  directions.  One 

can  write  the  local  force  equilibrium  equations  2.7  in  vector  form  as 

^  +  MTn  +  m#c+fc  =  0  (2.12) 

where  T,  and  G  refer  to  local  and  global  components.  Multiplying  by  [T)r  and  noting 
the  basic  kinematic  relations  given  by  equation  2.2  gives, 

[r\T^  +  ^f  n  + fa +  mTn=0  (2.13) 

and  upon  rearranging, 

dFr 

—  0  (2-14) 

where  one  has 

Pgt  =  PG  +  [T}TpL 

Fg  =  [T}t  Fl  (2.15) 

Fl  =  [T)Fg 
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In  scalar  form,  equation  2.14  becomes, 


dF, 

da 

dF\ 

da 

dj\ 

da 


=  ~PzT 


=  ~PvT 


=  ~PzT 


(2.16) 


Similarly,  one  can  write  the  local  moment  equilibrium  equations  2.8  in  vector  form 
as, 


dMi 

da 


0 

+  [/c]r  Ml  +  [T]  me  +  mi  +  ^  —  Fj 


=  0 


(2.17) 


Applying  the  same  transformations  as  for  the  force  equilibrium  equations  results  in, 


OMg 

da 


0 


-r  rhcT  +  [Tjr  s  — F 3 


=  0 


(2.18) 


where  one  has  defined, 


mCr  =  mc  +  [T]T  rhL 
Mg  =  [T]r  Ml 
Ml  =  [ T)Mg 


(2.19) 


In  scalar  form,  equation  2.18  becomes 

dMr 


da 

dMv 

da 

dM, 

da 


+  mzT  —  Tj\  Fi  +  T^\  F2  =  0 

+  rnvx  —  T22  ^3  +  ^32  F2  =  0 

+  mtj  —  T23  F3  4-  T33  F2  —  0 


(2.20) 


The  local  force  components  are  related  to  the  global  components  from  equations  2.15 


as 
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(2.21) 


Fi  =  Tu  Fz  4-  Tn  Fy  4-  T13  Fx 

F2  =  T21  Fx  4-  T22  Fv  4-  T23  Fx 

F3  =  T31  Fx  +  T32  Fy  +  T33  Fx 

One  places  the  above  into  equations  2.20  and  simplifies  by  using  the  following  relations 
obtained  from  noting  that  [T]_1  =  [T]T  and  applying  Cramer’s  rule  with  |T|  =  1. 


_  HP  /71  /Ti  rri 

11  —  -122  -*33  —  *23  -*32 

T\  2  =  T23  T31  —  T21  X33  (2.22) 

_  rri  /ti  rra  rr> 

13  —  12\  -*32  ~  -*22  ^31 


This  will  result  in  the  three  scalar  equations, 

dMt 


da 

dMv 

da 

dM x 


4-  mxr  —  T13  Fv  4-  T12  Ft  =0 
+  myj  4-  T,3  Fx  —  Tu  Fx  =0 
+  mzT  —  T12  Fx  4-  Tu  Fu  =0 


(2.23) 


Taking  the  derivatives  of  the  last  two  equations  and  introducing  the  force  equilibrium 
equations  2.16  gives 


d2Mv  dmvT  ,^('rp\,rp  p  ®T\\  _ 

~ds'  +  -dr  +  di[  "Fi)  +  T"T‘T~F-~d7  ~  0 

d2M,  drn.r  S  „  _  dTu 

~as~ +  ~ar  ~  J,(Tu  F,)  ~  Tu  P‘T  +  F“  ~ar  - 


(2.24) 


In  addition  to  these,  it  is  convenient  to  keep  the  local  moment  equilibrium  in  the  £ 
direction, 


dM  1 

da 


-  KC  M2  4 -  KnM3  +  mu  =  0 


(2.25) 


The  above  moment  equations  together  with  the  three  global  force  equations  2.16  are 
the  equivalent  of  equations  (71  b,  c)  (74),  and  (69  a,  b,  c)  of  Hodges  and  Dowell,  Ref. 
5.  No  approximations  have  been  made  as  yet  in  equations  2.16,  2.24,  2.25. 
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By  differentiating  equations  2.27,  then  solving  for  /?'  and  xp'  keeping  terms  only  to 
second  order,  one  obtains  the  same  expressions  as  would  have  been  obtained  by 
simply  differentiating  equations  2.28  directly.  Finally,  substituting  the  f3'  and  xp'  into 
the  three  curvature  strains  kj,  kp,  given  by  equations  2.4  and  keeping  terms  to 
second  order,  results  in, 
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k€  ~  8'  +  w'  v" 

Kn  —  v"  sin  8  -  w"  cos  8 

~  v"  cos  8  +  w"  Bin  6  (2.29) 

The  three  curvature  strains  are  now  expressed,  to  second  order,  in  terms  of  global 
deflections  v,  w  and  Euler  rotation  angle  8.  Often,  it  is  more  convenient  to  express 
the  twisting  behavior  of  the  blade  in  terms  of  a  total  twist  angle  <f>  which  is  defined 
as, 

<p  =  f  Kfda  =  8+  [  w' v" da  (2.30) 

Jo  Jo 

In  this  case,  the  curvature  strain  k ^  and  the  Euler  angle  8  are  replaced  in  equa¬ 
tions  2.29  by, 


*«  = 
8  = 


(2.31) 


Since  the  correction  to  the  Euler  angle  is  a  small  nonlinear  term,  it  is  often  neglected 
and  the  relation  8  ~  <f>  is  used. 

The  second  order  approximations  to  the  Euler  angles  as  given  by  equations  2.27 
are  also  used  for  the  general  transformation  matrix  [T].  Placing  these  trigonometric 
relations  into  the  basic  transformation  matrix  [T],  equation  2.1  gives  to  second  order, 


\T\  ~ 


\  1  -  vn!2  -  u/V2  v'  w' 

—  (v'  cos  8  -I-  w'  sin  8)  cos  8  (1  —  f'2/2)  sin  8  (1  -  uin/2) 

(a/ sin  6  —  w1  cos  8)  -  sin  8  (1  —  v^/2)  cos8(l-wn/2) 


(2.32) 
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The  third  step  in  the  reduction  process  is  to  relate  the  moment  resultants  to  the 
curvature  strains,  and  then  to  the  coordinates  v,  w,  9.  Using  the  generalized  linear 
stress-strain  relations  given  in  equation  2.11  and  introducing  the  strain-displacement 
relations  of  equations  2.29,  one  may  write, 


Aft 

—  E44  Kf 

~  GJ{9'  +  w'v ") 

m2 

=  £55  Hr, 

~  Elr,  (t/'sin  9  —  w" 

cos  9) 

m3 

—  Ees 

~  EIc{v"  cos  9  +  w" 

sin  9) 

The  above  are  for  a  blade  pricipal  axis  system  located  alon»  the  elastic  axis,  where 
there  is  no  coupling  between  the  (,  t],  and  (  axes.  For  non-principal  axes,  there 
may  be  additional  couplings  between  tj  and  £  and  for  non-elastic  axis,  such  as  in 
composite  blades,  there  may  be  additional  couplings  between  the  £  and  rj  and  £  and 
(  curvatures.  For  use  in  the  equilibrium  equations  2.24,  it  is  also  required  to  express 
the  moments  in  global  x,  y,  z  directions  in  addition  to  the  local  £,  77,  £  directions 
given  by  equations  2.33.  From  equations  2.19,  one  has 

Mx  =  TnM\  +  T21M2  +  T31M3 

Mv  =  T[ 2 A-f  1  4-  T22A/2  4-  T32A/3  (2.34) 

Mt  =  Ti3iM\  4-  T23A/2  4-  T33A/3 

This  gives,  to  second  order, 

Mv  =  GJ6'v’-{EIt  sin20  +  £/„  cos20)u/' 

-{EI{  —  EIn)  cos  9  sin  9v" 

Mt  =  GJ  9'  w'  4-  (£7C  cos2  B  4-  EI„  sin2  6)  v"  (2.35) 

+(EI(  -  EIn)  cos  9  sin  9  w" 

Mx  is  not  given  above,  since  in  the  present  formulation,  the  local  moment  M\  is  used 
rather  then  the  global  moment  Mx. 
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Finally,  to  complete  the  reduction  process,  one  places  the  moments  equations  2.35, 
2.33  and  curvature  strains  equation  2.29  into  the  equilibrium  equations  2.24,  2.25  to 
obtain, 

[GJe'v'  -  (£/c  sin2  8  +  EI„  cos 78)w"  -  ( EIC  -  EI„)  cos  8  sin  6v”}" 

+(u/  Ft)'  +  (1  -  v'7 /2  -  w'7 /2)ptr  -  Fx  (v  v”  +  w’  w")  +  m'yT  =  0 

[GJ  6'  w'  +  ( EIC  cos2  8  +  El,,  sin2  8)  v"  +  (£/c  -  £/„)  cos  8  sin  8  w"\" 

~{v'  Fzy  +  (1  -  vr*/2  -  w'7 /2)pvr  +  Fu  (v'v"  +  w'  w")  +  m'lT  =  0  (2.36) 

[GJ  ( 8 '  +  w'  v")}'  -  {EIC  ~  i(^"2  -  U"2)  cos  8  sin  8  +  v"w"  cos  28} 

+m it  -  0 

The  force  loadings  Fz ,  Fv,  Ft  in  the  above  are  found  from  integrating  the  global  force 
equations  2.16.  For  free  vibrations,  the  inertial  loe-uus  p,,  pv,  px  and  m i  are  given 
by  equations  2.9  and  2.10. 

Although  the  above  equation^  have  been  reduced  formally  to  second  order,  some 
further  simplicaticas  are  still  made  to  reduce  them  to  a  simpler  form.  First,  as 
mentioned  in  Ref.  4,  by  integrating  the  third  equation,  then  multiplying  it  by  v ',  then 
subtracting  it  from  the  first  equation,  one  can  eliminate  the  GJ  8'  v'  term,  introducing 
only  new  third  order  terms  from  the  third  equation.  Hence,  to  second  order,  the 
GJ  8'  v'  term  can  be  neglected.  Similarly  for  the  GJ  8'  w'  term  in  the  second  equation. 
Next,  the  v n  and  w n  terms  can  be  neglected  compared  to  unity  for  moderate  deflection 
slopes.  This  would  also  eliminate  the  Fu  and  Fx  terms  since  they  were  multiplied  by 
and  T\\  is  now  set  equal  to  unity  as  seen  in  equation  2.32.  Along  the  same  lines, 
all  derivatives 

d_ 

da 

in  these  equations  can  be  replaced  by 

d_ 

dx 
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since  from  the  kinematic  relations,  equations  2.6, 

d  ddx  n  n  d  d 

!2~V  l2^di~dl 


(2.37) 


Also,  it  is  convenient  to  introduce  the  total  twist  variable  ^  as  defined  by  equation  2.30 
rather  than  deal  with  the  Euler  angle  9.  With  these  simplifications,  the  previous 
eqautions  can  be  rewritten  as, 


w  :  sin2  9  +  EIn  cos2  9)  w"  -f  ( EIq  —  EIn)  cos  9  sin  9  v"}" 


~(w'  Fz)’  =  pzT  +  m' 


vT 


v  :  [( EI{ i  cos2  9  +  EI^  sin2  9)  v"  +  ( El (  —  EIn)  cos  9  sin  9  w")" 

-WF,y  =  PvT-m'vT  (2-38) 

O  :  -( GJ  <t>')'  4-  ( EI^  —  EI^)  f(u/'2  —  v”2)  cos  9  rin  9  4-  v"w"  cos  29} 

=  rriiT 


where  one  has, 

9  ~  <f>  —  [  w1  v"  dx 
Jo 

Fx  ~  +  J  PzT  dx  (2.39) 

e  -  u'  +  v'2/2  +  wn/2  =  0 


Equations  2.38  are  effectively  the  nonlinear  moderate  deflection  equations  presented 
by  Hodges  and  Dowell,  Ref.  4,  Boyd,  Ref.  5,  and  others.  They  have  been  shown 
to  arise  from  a  straightforward  reduction  of  the  general  nonlinear,  large  deflection 
equations  given  by  Minguet  and  Dugundji,  Ref.  1,  and  presented  here  in  section  2.1. 
Often,  the  relation  9  ~  0  is  used  in  place  of  the  more  accurate  relation  given  by 
equations  2.39.  The  e  =  0  relation  of  equations  2.39  represents  an  effective  no  stretch 
condition  and  is  used  to  determine  the  axial  deflection  u  since  v  and  w  deflections 
have  been  determined.  For  vibration  problems  the  inertia  loadings  are  given  by 
equations  2.9  and  2.10  with  -Ip9  replaced  by  —  Ip<f>. 
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One  last  item  of  reductions  of  these  equations  comes  about  by  eliminating  the 
trigonometric  functions  in  them.  For  a  flat  blade  without  built-in  twist,  8t  =  0,  the 
trigonometric  functions  can  be  expanded  to  second  order  as, 

sin  6  ~  8 

cos0  ~  l-02/2  (2.40) 

Placing  these  into  the  previous  equations  2.38  gives  the  more  useful  form, 

w  :  [( EI„w "  +  (EIC  -  EI„)  (t ,"8  +  w"  82)\" 

~(w'  Fx)'  =  pzT  +  m'r 

v  :  [ EIcv "  4 -(£/<-  El^iu,"  8  -  v"  82)}" 

~(v' Fx)' =  pvT  -  m'xT  (2.41) 

< t>  :  ~{GJ  <t>' )'  +  {EIC-  Elr,)  [(u;"2  -  v"2)  8  +  v"  w"] 

=  m\T 

This  form  shows  more  clearly  the  type  of  nonlinear  couplings  involved  between  the 
w ,  v ,  and  0  motions.  These  nonlinear  couplings  depend  on  the  difference  in  bending 
stiffness,  ( El <  —  EIn),  and  would  give  rise  to  linear  couplings  by  the  presence  of  an 
initial  static  deflection  in  w  and  v.  Similar  equations  can  be  obtained  for  blades  with 
an  initial  twist  8t,  by  replacing  equations  2.40  with, 

sin(0t  +8)  ~  sin  8t  -1-  8  cos  8t  —  \f?2/2')sin  8t 

cos(8t  +  8 )  ~  cos 8t  —  8  sin#<  —  [82/2\cosdt  (2.42) 

Although  the  moderate  deflection  equations  2.41  lend  themselves  well  to  Galerkin 

solution,  one  should  be  careful  to  use  a  sufficient  number  of  modes  to  capture  the 

nonliner  effects  when  static  deflections  are  present.  They  can  always  be  checked 
against  the  general  solution  of  the  twelve  nonlinear  differential  equations  presented 
by  Minguet  and  Dugundji,  Ref.  1. 
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Before  leaving  this  section,  it  might  be  interesting  to  note  that  the  moderate 
deflection  equations  can  also  be  derived  from  an  energy  formulation  by  minimizing 
the  total  potential  energy  II  of  the  functional, 

II  =  -  /  EIn  (w"  cos  9  -  v"  sin  9)2  dx 
2  Jo 


+  -  f  EI(  (w"  sin  9  +  v"  cos  6)2  dx 
2  Jo 

+  1  fL  GJ  (6'  +  w'v")2dx 
2  Jo 

+  ]-  [  Fx(w'2  +  v'2)dx 
2  Jo 

-  /  ( PUT  V  +  PiT  W  -  muT  v'  +  ™zT  v'  +  m\T  <f>)  dx 

Jo 


(2.43) 


A  simple  application  of  variational  methods  will  lead  to  the  moderate  deflection  equa¬ 
tions  given  by  equations  2.38  and  2.41. 


2.3  Modeling  of  Large  Amplitude  Motion 

In  Ref.  1,  the  basic  equations  given  in  section  2.1  were  linearized  around  a  given 
static  position  to  yield  a  small,  perturbed  free  vibration.  An  appropriate  eigenvalue 
problem  was  then  solved  to  find  the  various  mode  shapes  and  their  associated  natural 
frequencies.  This  eigenanalysis  is  not  useful  for  large  amplitude  motion  because  once 
structural  nonlinearities  are  present  in  both  static  and  dynamic  sense,  the  natural 
frequency  of  a  particular  mode  becomes  a  function  of  amplitude  of  that  mode.  Fur¬ 
thermore,  it  is  also  expected  that  certain  amount  of  couplings  exist  between  the  static 
and  dynamic  components  in  the  various  variables.  Thus  two  basic  characteristics  that 
distinguish  the  nonlinear,  large  amplitude  vibration  from  the  linear,  small  vibration 
can  be  summarized  as  follows. 

(1)  The  natural  frequency  of  a  particular  mode  changes  as  its  amplitude  increases. 

(2)  The  static  mean  position  of  the  beam  can  also  change  as  a  function  of  amplitude. 
Two  popular  methods  for  the  solution  of  general  nonlinear  dynamic  problems  are 

direct  numerical  time  integration  of  the  basic  equations,  and  the  harmonic  balance 
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method.  The  former  method  will  give  the  exact  solution  which  shows  the  effects 
of  all  possible  harmonics,  while  the  latter  method  will  yield  a  solution  with  only 
few  harmonics.  The  direct  time  integration  requires  a  set  of  governing  equations 
that  contain  only  time  t  as  independent  variable.  By  performing  appropriate  modal 
analysis,  one  has  to  reduce  the  equations  of  motion  into  a  modal  form,  expressing 
them  as  functions  of  generalized  coordinates.  Usually,  a  large  amount  of  computing 
time  is  used  until  the  solutions  reach  their  final  steady  states. 

In  the  present  analysis,  the  harmonic  balance  method  is  used  because  we  do  not 
want  to  begin  with  a  set  of  approximate  modal  equations  which  are  based  on  an 
ordering  scheme,  but  rather  use  the  large  deflection  equations  of  section  2.1.  These 
twelve  differential  equations  contain  all  the  twelve  variables,  i.  e.  three  Euler  angles, 
three  force  resultants  and  three  moment  resultants,  in  addition  to  the  usual  three 
displacements  x,  y,  z  as  their  independent  variables.  In  such  a  situation,  it  is  more 
insightful  to  assume  the  time  dependency  of  the  solution  in  the  the  form  sinu;f,  and 
use  numerical  integration  in  space  rather  than  in  time.  In  doing  so  one  loses,  of 
course,  the  effects  of  higher  harmonics,  but  the  key  arguement  is  that  in  most  of 
the  nonlinear  analysis,  amplitudes  associated  with  the  first  harmonic  take  the  largest 
quotient, and  therefore  are  most  critical  in  determining  its  response  or  stability. 

Thus  for  the  purpose  of  present  analysis,  all  quantities  are  assumed  to  be  of  the 
following  form 

X(w,  f )  =  Xo{u>)  +  X,{v)  sin  u it 

where  Xq,X,  represent  the  static  part  and  the  associated  amplitude  (not  a  small 
quantity)  around  that  static  part,  respectively.  The  fact  that  X,  is  not  a  small 
quantity  is  reflected  in  the  frequency  dependency  of  both  X0  and  X,.  Hence,  unlike 
small  vibration  problem,  there  exists  one-to-one  correspondence  between  amplitude 
and  frequency. 

The  analytic  modeling  consists  of  substituting  the  above  expression  for  each 
variable  into  the  twelve  governing  equations.  As  a  result  of  multiplications  involving 
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sin  u/t,  this  will  produce  many  higher  order  terms  containing  higher  harmonics  such 
as  sin  2 ut,  sin  3a it.  For  details  of  how  these  multiplications  are  performed  and  what 
the  resulting  coefficients  are,  see  the  Appendices.  A  harmonic  balance  method  is 
then  employed  to  retain  only  two  kinds  of  terms;  the  ones  that  are  constants  and 
the  ones  that  are  coeffcients  of  sinwt.  All  the  higher  harmonic  terms  are  left  out. 
Some  of  the  remaining  terms  will  contain  higher  order  of  magnitude  terms,  for  ex¬ 
ample,  sin4  ut  produces  the  constant  3/8  even  after  neglecting  its  higher  harmonic 
components  cos  2ujt  and  cos4a;f.  It  is  clear  that  keeping  all  these  higher  order  of 
magnitude  terms  will  make  the  equations  extremely  long  and  unwieldy.  Hence,  an 
ordering  scheme  that  keeps  magnitudes  of  up  to  third  order  is  employed  to  maintain 
a  consistent  level  of  nonlinearities  in  all  of  the  equations.  See  the  Appendices.  It  is 
emphasized  that  this  ordering  scheme  does  not  mean 

Q 2 

cos  9  ~  1  — —  +  H.  0.  T. 

2 

but  rather 


cos#  ~  cos  #0  ~  sin #o  A#  —  -  cos  #o  ( A#)2 

+  ^  sin  #0  (A#)3  +  H.  0.  T. 

6 

where  9  =  #0  4-  A 9,  and  the  90  and  A 9  —  93sinu>t  represent  the  static  and  dynamic 
component*  of  9.  So  the  complete  nonlinearity  in  the  large  rotations  is  still  kept  in  a 
static  sense,  but  a*  a  strategy,  terms  only  up  to  third  order  are  kept  in  the  dynamic 
counterpart*. 

One  point  is  noted  here;  applying  the  harmonic  balance  followed  by  the  approxi¬ 
mating  schemes  will  not  render  the  final  twenty  four  equations  completely  compatible 
with  each  other.  More  specifically,  these  coupled  equations  would  not  satisfy  equilib¬ 
rium,  geometric  compatibilities,  and  stress-strain  relations  perfectly  as  their  original 
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twelve  versions  would.  Therefore,  one  should  expect  deterioration  in  the  degree  of 
compatibility  as  amplitudes  increase.  Normally  this  would  mean  loss  of  accuracy  in 
the  solutions,  or  in  the  worst  case,  even  the  loss  of  convergence.  However,  as  shown 
later  in  this  report,  this  does  not  impose  serious  computational  limits  in  most  of 
practical  range  of  amplitudes. 
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Chapter  3 

Method  of  Solution 


Having  obtained  all  the  necessary  formulations  for  the  large  amplitude,  nonlinear  free 
vibration  model,  one  can  express  the  equations  of  section  2.1  in  vector  form 

~  =  *,(*„, X„uO  (3.1) 

(12x1)  (12x1) 

and 

^  =  g.(X0,X.,u)  (3.2) 

(12x1)  (12x1) 

where 

*0  ~  [  Fio  Fio  Fx  A/ io  A/2o  A/30  Xo  yo  zq  60  /?o  V'o  ] r 

X,  =  [F\,  F2,  F3,  M\,  M2.  M3,x,y,  z,9,/3,ip,}T 

The  two  vector  function  arrays  go  and  g,  contain  many  product  terms  involving 
multiplications  of  two,  or  three  harmonic  quantities.  They,  of  course,  originate  from 
the  twelve  basic  equations  that  are  presented  in  section  2.1.  Multiplications  of 
harmonics  and  calculations  of  the  coefficients  of  the  resulting  new  harmonics  can  be 
easily  implemented  according  to  the  formulae  in  the  Appendices. 

To  solve  this  system,  all  of  the  twenty  four  equations  (now  twelve  for  the  static 
part,  twelve  for  the  dynamic  part)  are  first  integrated  from  the  tip  to  the  root  of  the 
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blade  once.  In  the  previous  TEL  AC  report  Ref.  1,  Minguet  used  a  finite- difference 
iteration  method  for  the  solution  of  static  deformation,  sweeping  from  the  tip  to 
the  root  and  vice  versa  a  few  times  until  all  the  residues  become  very  small.  When 
applying  this  scheme  to  the  solution  of  mode  shapes  and  their  frequencies,  one  has 
to  be  cautious  because  this  finite-difference  iteration  will  usually  converge  to  the  first 
mode  only.  To  obtain  higher  modes,  one  must  consider  other  integration  techniques 
which  do  not  sweep  back  and  forth  along  the  span  but  are  more  appropriate  for 
boundary  value  type  problems.  Among  such,  Runge-Kutta  integration  is  frequently 
used  and  very  effective.  Currently  fourth  order  Runge-Kutta  algorithm  is  used. 

In  the  early  step  of  numerical  integration,  one  has  to  guess  boundary  values  of 
displacements  and  rotations  at  the  tip  as  well  as  the  frequency  that  will  make,  for 
a  given  mode  shape,  all  the  displacements  and  rotations  at  the  root  as  close  to  the 
prescribed  values  as  possible.  For  instance,  a  linear  solution  by  Minguet  can  provide 
such  a  good  guess  for  tip  values  Xt.  The  functional  relationships  between  these  two 
sets  of  boundary  values  at  the  root  and  at  the  tip  can  be  written  as 

Xr  =  f(X„u)  (3.3) 

(12x1)  (12x1) 

where 

Xt  =  [xoz.yoy,  z0z,  8o9t{3Q0,  V'o^,]7' 

at  the  tip,  and 

Xr  =  [0  0  0  0  0  0  0O  0  /30  0  V’o  0]T 

at  the  root. 

Here  8o,0u,^o  are  prescribed  values  at  the  root  (they  are  zero  for  flat  cantilever 
blades).  Since  the  initial  guess  for  the  twelve  components  of  Xt  can  not  be  perfect, 
there  will  be  nonzero  residues  R  by  the  time  the  integration  reaches  the  root.  A 
Newton-Raphson  type  algorithm  can  then  be  used  to  produce  a  better  set  of  boundary 
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values  based  on  the  current  values.  This  will  produce  a  series  of  the  following  set  of 
boundary  values. 


Xtn+1  =  X t"  -  J(Xtn,vn)-lRn 

where 


(3.4) 


ft"  =  f(X?,u>~)  -  X, 

and 


J  :  (12  x  12)  Jacobian  matrix 

Here  the  superscript  n  refers  to  the  n  th  iterative  values,  and  Xr  refers  to  the  desired 
values  at  the  root.  The  n  th  boundary  values  X"  at  the  tip  will  eventually  march  to 
the  true  solution,  provided  it  exists.  Currently  two  algorithms  called  F.  D.  G.  (finite 
difference  Gauss’  method)  and  F.  D.  L.  M.  (finite  difference  Levenberg-Marquardt 
method)  (Ref.  7),  respectively  are  used.  The  former  is  simply  a  numerical  version 
of  Newton- Raphson  method,  and  in  the  latter  case,  an  efficient  relaxation  scheme  is 
added. 

It  is  noted  that  whatever  algorithm  is  used,  it  must  take  iterations  on  the  frequency 
as  well  as  the  boundary  values,  since  it  is  not  known  in  advance  at  which  frequency 
a  mode  will  happen  for  a  given  amplitude  level.  Therefore,  one  of  the  six  boundary 
amplitudes  at  the  tip  x,,y„  is  replaced  by  the  frequency  u>,  and  the 

replaced  displacement  is  fixed  throughout  iterations.  Which  one  has  to  be  fixed 
depends  on  which  mode  is  being  sought.  For  instance,  if  bending  modes  are  of  concern 
it  will  be  z,;  if  it  is  torsional  modes  then  9,  is  fixed.  The  iteration  will  march  until 
the  boundary  conditions  at  the  root  are  met,  i.  e.  ,  the  residues  Rn  are  zeros  or  at 
least  less  than  some  preset  e  where  e<  1. 

As  a  final  notion,  the  above  solution  procedure,  when  applied  to  linear  problems, 
is  similar  to  the  so  called  transfer  matrix  technique  used  to  obtain  helicopter  blade 
vibration  modes  by  Isakson  and  Eisley  in  Ref.  8. 
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Chapter  4 

Results  and  Discussion 


The  prescribed  algorithms  have  been  used  to  investigate  the  first  and  second  bending 
modes,  first  fore-and-aft  modes,  and  first  torsion  modes  of  cantilevered  blades  with 
the  lay-ups  [0/90]3,,  [45/0],  of  graphite/epoxy  for  various  tip  deflections.  These  modes 
were  chosen  because  tv  -•>  have  the  lowest  natural  frequencies  and  hence  should  be 
easily  converged.  aermore,  they  pose  much  importance  from  a  aeroelastic  point 
of  view.  The  c  ^figurations  of  the  blades  investigated  are  those  used  by  Minguet  in 
Ref.  1  (560  mm  long,  30  mm  wide).  Beam  material  properties  of  these  lay-ups  are 
listed  on  Table  1.  To  see  how  these  coefficients  are  calculated,  refer  to  section  2.6  of 
Ref.  1. 

The  static  deflections  were  varied  by  imposing  and  adjusting  uniform  gravity  level 
throughout  the  blade.  As  stated  earlier,  one  of  the  six  boundary  amplitudes  at 
the  tip  was  replaced  by  o>,  and  the  replaced  amplitude  was  fixed  throughout  the 
iterations.  The  z,,y,,  and  9,  were  fixed  for  bending,  fore-and-aft,  and  torsional  modes, 
respectively.  Also,  9o,/3o,iI>q,  at  the  root  were  all  set  equal  to  zero  since  the  blade 
is  a  flat  cantilever  beam.  A  total  of  16  node  points  were  used  along  the  blades. 
Note  that  the  same  number  of  nodes  was  also  used  in  Ref.  1.  All  of  the  cases  were 
guided  by  the  linear  results  by  Minguet.  That  is,  the  linear  mode  shapes  and  their 
natural  frequencies  provide  reasonable  trial  values  which,  after  a  few  iterations,  would 
lead  to  nontrivial  solutions.  All  the  runs  were  made  on  a  DEC  Microvax  computer 
with  typical  number  of  iterations  from  5  to  10  for  convergence.  Each  iteration  took 
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approximately  between  15  to  30  seconds  of  CPU  time,  longer  times  being  required 
for  cases  with  strong  structural  couplings.  Very  often,  it  was  necessary  to  use  under- 
relaxation  to  lead  iterations  smoothly  to  the  final  solution  without  causing  divergence 
or  any  sudden  jump  into  another  nonlinear  solution  (In  fact,  both  F.  D.  G.  and  F. 
D.  L.  M.  algorithms  assume  use  of  certain  under-relaxations).  Each  analysis  was 
continued  until  the  amplitude  could  not  be  further  increased.  At  this  point,  the 
Jacobian  matrix  became  almost  singular  and  the  solution  did  not  converge. 

Before  illustrating  the  results  in  detail  it  is  worthwhile  to  mention  that  in  linear 
problems  where  perturbations  are  very  small,  the  present  analysis  would  be  slightly 
superior  to  Ref.  1.  The  present  analysis  is  based  on  a  continuous  model  while  Ref.  1 
is  based  on  a  lumped,  finite  difference  model. 

The  first  example  is  that  of  the  [Q/90j3,  specimen.  Figure  3  through  figure  26 
show  mode  shapes  at  two  amplitude  levels  under  three  different  static  tip  deflections. 
It  is  seen  that  for  most  of  the  amplitude  range,  the  nonlinear  modes  remain  almost 
the  same  as  linear  modes  in  their  shapes  even  though  their  frequencies  change.  Next, 
figure  27  through  29  show  change  of  natural  frequencies  as  functions  of  amplitudes 
2,,y,,  and  9,  at  the  tip.  Also  Figure  30  and  31  represent  the  variations  of  the  cen- 
tershifts  z0  at  the  tip  of  various  modes  as  functions  of  the  tip  amplitudes.  From  the 
figures  the  following  two  observations  can  be  made. 

(1)  Increasing  amplitude  level  has  slight  stiffening  effects  in  IB,  2B  (or  any  bending 
modes,  presumably)  whereas  it  has  significant  softening  effects  in  IF,  IT  modes, 
particularly  for  moderate  range  of  static  tip  deflections.  As  a  result,  the  natural 
frequencies  of  bending  modes  rise  slightly  with  amplitude  level  while  those  of  IF  and 
IT  modes  always  drop. 

(2)  The  above  frequency  changes  are  accompanied  by  centershift  changes.  Increasing 
amplitude  levels  has  slight  effects  on  the  centershifts  of  bending  modes  except  for 
the  2B  mode,  whereas  it  has  significant  centershift  increase  for  the  IF  mode  and  a 
centershift  decrease  for  the  IT  mode,  particularly  for  moderate  static  tip  deflections. 
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The  behavior  of  these  centershift  changes  seem  relavant  to  the  linear  findings  in  Ref. 
1  (see  figure  33). 

Figure  32  presents  the  effects  of  second  harmonics  on  the  natural  frequencies  of 
IF  modes.  It  was  found  that  including  terms  involving  the  second  harmonic  cos  2 u>t 
in  two-dimensional  sense  was  enough  to  capture  the  missing  second  harmonics  in  IF 
modes.  In  other  words,  only  F\,Fz,M 2,  and  0,x,z  were  expressed  in  the  form 

X(u,  t )  =  Xo(w)  -I-  X.{u j)  sin  u>t  +  Xic(v)  cos  2 wt 

with  all  other  variables  containing  only  the  first  harmonics  as  before.  This  was  done 
based  on  the  intuition  that  second  harmonics  will  mostly  appear  in  z„  z,  and  their 
motion  should  be  initially  90  degrees  out  of  phase  with  the  rest  of  amplitudes.  Then 
a  new  set  of  formulae  that  performs  multiplications  of  harmonics  was  implemented 
in  the  computer  program.  These  are  different  from  the  previous  ones  in  the  appendix 
in  that  they  now  have  to  deal  with  cos  2 uit  as  well.  The  resulting  Jacobian  is  then 
(15  x  15)  instead  of  (18  x  18)  which  would  result  if  cos  2 ut  were  introduced  in  ail 
of  the  variables.  As  can  be  seen  from  the  plots,  IF  modes  exhibit  significant  second 
harmonic  contents  in  z,  motion  for  moderate  range  of  static  tip  deflections  (roughly, 
from  20  mm  to  80  mm.)  as  amplitude  is  increased.  On  the  other  hand,  at  either  zero 
or  very  large  tip  deflection  the  second  harmonics  are  almost  unrecognizable.  In  fact, 
z$  has  no  first  harmonic  content  in  IF  modes.  An  effort  was  also  made  to  seek  for 
any  second  harmonics  in  IT,  IB  and  2B  bending  modes,  but  they  have  been  found 
very  weak  and  are  not  presented  here. 

Next  example  is  that  of  [45/0],  which,  unlike  the  previous  case,  exhibit  bending- 
torsion  coupling.  Due  to  the  structural  coupling,  computer  time  was  increased  and 
the  convergence  became  more  sensitive.  This  resulted  in  earlier  breakdown  of  nonsin¬ 
gularity  of  Jacobian  matrix  which  in  turn  caused  shorter  range  of  solutions  available 
as  functions  of  amplitudes.  Figure  34  through  57  show  mode  shape  changes  at  two 
amplitude  levels  under  three  different  tip  deflections.  Once  again,  the  mode  shapes 
do  not  change  significantly  from  the  linear  modes.  Figure  58  through  62  show  the 
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frequency  and  centershift  changes  as  the  amplitudes  of  various  modes  increase.  The 
two  former  observations  (1)  and  (2)  can  also  be  made  in  these  figures;  a  similar  anal¬ 
ogy  about  the  relationship  between  frequency  and  centershift  changes  can  be  also 
made.  The  effects  of  second  harmonics  on  the  natural  frequencies  of  IF  modes  is 
shown  in  figure  63.  Unlike  the  previous  case  of  [Q/90]3f,  the  presence  of  second  har¬ 
monics  is  relatively  weak.  In  particular,  due  to  the  existing  bending-torsion  coupling, 
the  static  tip  deflection  will  not  lie  on  the  z  axis,  and  the  IF  motion  is  not  symmetric 
about  the  z  axis  even  though  the  root  angles  here  are  again  zeros. 

Finally,  it  is  interesting  to  consider  what  makes  the  Jacobian  matrix  singular  at  a 
certain  point  along  the  way  of  increasing  amplitude.  Except  for  the  cases  of  IB,  there 
seem  to  be  certain  limits  on  the  largest  amplitudes  that  can  be  solved  by  the  current 
algorithms.  These  limits  were  even  more  severe  if  second  harmonics  were  included. 
In  section  2.3,  it  was  suggested  that  one  should  expect  deterioration  in  the  degree 
of  compatibility  as  amplitudes  increase.  This  could  be  one  possibility.  Apart  from 
that,  other  factors  may  attribute  to  the  singularity  of  solution;  the  round-off  errors 
associated  with  the  large  size  of  Jacobian  matrix,  and  the  interaction  of  several  modes 
as  amplitudes  increase,  with  possible  resulting  chaotic  vibration. 
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Chapter  5 
Conclusion 


Throughout  the  research  period,  it  has  been  demonstrated  that  the  nonlinear  analysis 
derived  from  work  by  Minguet  in  Ref.  1,  and  iteration  methods  based  on  harmonic 
balance  and  numerical  integration  of  the  basic  equations  seems  efficient  for  large 
amplitude  vibration  problems  of  composite  rotor  blades.  These  include  the  nonlinear 
free  vibration  problem  which  is  presented  here,  and  nonlinear  limit  cycle  problems 
with  dynamic  stall  in  hover  and  possibly  in  forward  flight,  which  will  be  investigated 
as  parts  of  future  work. 

For  the  free  vibration  part,  it  has  been  shown  that  both  large  static  deflection  and 
large  amplitude  can  affect  significantly  the  fore-and-aft  modes  and  torsion  modes,  but 
not  much  the  bending  modes.  More  specific  conclusions  are  as  follows. 

(1)  Increasing  amplitude  level  has  slight  stiffening  effects  in  bending  modes  whereas 
it  has  significant  softening  effects  in  IF,  IT  modes,  particularly  for  moderate  range 
of  static  tip  deflections.  As  a  result,  the  natural  frequencies  of  bending  modes  rise 
slightly  while  those  of  IF  and  IT  modes  always  drop. 

(2)  Increasing  amplitude  level  of  a  particular  mode  also  results  in  centershift  changes 
that  are  small  for  the  bending  modes  but  significant  for  the  IF  and  IT  modes, 
particularly  for  moderate  static  tip  deflections.  The  IF  centershift  seems  to  increase 
considerably  with  amplitude  level.  The  behavior  of  these  centershift  changes  seem  to 
stem  from  the  linear  findings  in  Ref.  1. 

(3)  The  flat  [90 /0] 3,  or  any  isotropic  blade  with  zero  root  angle  has  significant  second 
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harmcaic  contents  in  the  IF  mode  for  moderate  static  tip  deflections.  These  appear 
mostly  in  the  bending  amplitude  z,.  If  the  root  angle  is  not  zero,  or  there  is  bending- 
torsion  coupling  however,  the  second  harmonics  are  not  as  strong. 

The  conclusions  made  during  the  research  period  should  give  insight  into  more 
general  nonlinear,  large  amplitude  analysis  such  as  proposed  by  Dugundji  (Ref.  9)  for 
future  investigation. 

Regarding  the  future  work  which  is  specified  in  Ref.  9,  nonlinear  limit  cycle  anal¬ 
ysis  of  composite  blades  in  the  presence  of  dynamic  stall  is  currently  being  pursued. 
The  structural  nonlinearities  are  well  represented  by  the  current  model,  and  for  the 
aerodynamic  nonlinearity,  the  ONERA  model  developed  by  Tran  and  Petot  in  Ref. 
10  is  used.  The  analysis  begins  with  a  simple  two-dimensional  motion  with  bending 
modes  only  (called  “bending  stall”),  and  later  will  go  into  three-dimensional  motion 
with  additional  torsion  and  fore-and-aft  modes  as  well  as  centrifugal  forces,  Coriolis 
acceleration,  and  coning  angles  also  present. 
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Appendix  A 

Calculation  of  Coefficients  of 
Harmonic  Quantities 


In  section  2.3  it  was  suggested  that  for  large  amplitude  motion,  every  variable  be 
expressed  as 


X(u >,  t )  =  Ao(w)  +  X,(u>)  sin  u>t 

where  Xq,X»  represent  the  static  and  dynamic  components  of  a  particular  variables. 
As  a  result,  all  the  quantities  in  the  original  twelve  governing  equations  will  take  the 
above  form  immediately.  Recall,  however,  that  many  of  the  terms  in  the  equations 
involve  trigonometric  functions  and  their  arguments  are  the  three  Euler  angles  v,/3,9. 
Then  it  is  clear  that  one  can  not  apply  harmonic  balance  method  with  the  Euler 
angles  expressed  as  above  and  themselves  inside  the  trigonometric  functions.  So,  it  is 
useful  to  rely  on  series  expansion  versions  of  these  trigonometric  functions.  In  order  to 
get  the  series  expressions,  let  z  represent  any  of  the  three  Euler  angles,  and  let  A” ( x ) 
be  any  trigonometric  function,  i.  e.  cos  z,  sinz,  tanz,  or  1/cosz.  Then  substituting 

z  =  Zo  +  x,  sin  u>t 

into  the  function  X  and  expanding  in  a  Taylor  series  about  xq  yields 
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X(x)  =  X(x0)  +  —  (xo)x,  smut 

CIS 

(PX  (P  X 

+  1/2!  — ^-(x0)  x]  sin2  Ult  +  1/3!  ~ -(x0)  x]  sin3u it  (A.l) 

ax *  aiJ 

=  Xq  +  X,  sin  u/t 

+Xj2  sin2  ujt  +  X,3  sin3  u/t  +  H.  0.  T. 

where 


Xo 

=  X(xo) 

(A. 2) 

dX ,  v 

X, 

(A. 3) 

<PX  , 

x,2 

s  I/*  ^ 

(A. 4) 

X,3 

s  i/®  £<*.)»; 

(A.5) 

Here  according  to  our  ordering  scheme  only  terms  up  to  third  order  are  kept  in 
the  expansion  (see  section  2.3).  Then,  when  applying  harmonic  balance  methods, 
the  sin2^  and  sin3  ut  can  be  expanded  into  constant  and  sinud  type  terms  after 
multiplication  with  other  harmonic  quantities,  as  shown  in  Appendices  B  and  C. 

In  the  current  analysis  four  different  trigonometric  functions  are  encountered.. 
They  are  cos  x,  sinx,  tan  x,  and  1/ cosx.  According  to  above  expansion  rules  then 
each  trigonometric  function  can  be  expressed,  up  to  third  order,  as 


cos  x 


cos  xo  —  (sin  xo)  x,  sin  wt  -  1/2  (cos  xo)  x2  sin2  u/t 
+  1/6  (sin  x0)x3  sin3wt 


(A.6) 


sinx  =  sin  xo  +  (cos  xo)  x,  sin  u/t  —  1/2  (sin  x0)  x2  sin2  u/t 
— 1/6  (cos  x0)  x3  sin3u/£ 


(A.7) 
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tan  x 


1  /  cos  x  — 


=  tan  xQ  +  (1/  cos2 x0)  x,  sin  ut  -I-  (tan  xq/  cos2  Xo)  x2  sin2  ujt 

+  1/3  ((2  tan*  xo  +  1/  cos2  x0)  /  cos2  x0)  13  sin3  u >t  (A. 8) 


1/  cos  x0  +  (tan  x0/  cos  x0)  x,  sinu;t  +  1/2  (1  /cos2x0  +  tan2  x0/  cos  x0) 

.  x2  sin2  uit  +  1/6  (5  tan  x0/  cos3  x0  +  tan3  x0/  cos  x0)  x3  sin3  ut  (A. 9) 
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Appendix  B 

Multiplication  of  Two  Harmonic 
Quantities 


In  the  previous  section  it  was  seen  that  any  harmonic  quantity  can  be  expressed,  up 
to  third  order,  as 

X  =  Xo  +  X,  sin  ut  +  X,2  sin2  ut  +  X,3  sin3  ut  ( B .  1 ) 

where  X0,X,,X,2,X,3  are  determined  by  the  formula  A. 2.  X{x)  could  be  either  a 
harmonic  variable  itself  (e.g.  F,\f,x,...  etc.  )  or  a  trigonometric  function.  If  it  is  a 
harmonic  variable  X,2,X,3  are  identically  zero.  Now  let’s  consider  a  product  of  two 
quantities,  X  and  Y  which  are  expressed  as  above.  It  can  be  shown  that 

X  Y  —  {Xo  +  X 4  sin  ut  +  Xt2  sin2  ut  -+-  A ,3  sin3  ut ) 

.  (Vo  +  Y,  sin  ut  4-  Y,2  sin2  ut  +  Y,3  sin3  ut) 

=  {XY)0  +  (XY),  sin  ut  +  {XY),2  sin2  +  {XY)s3  sin3  ut  (B.2) 

where 


{XY)  0  =  A0y0 
{XY).  =  XoY.  +  X.Yo 
{XY),  2  =  X0Y,2  +  X.YS  +  X,2Vo 
{XY),3  =  A0y.3  +  X.Y.2  +  X,2Y,  +  X.3Yo 
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When  applying  harmonic  balance  method  only  the  static  and  the  first  harmonic  terms 
are  retained.  For  this  purpose  note  that 

sin2u/i  =  1/2  —  1/2  cos 2u>t  (B.3) 

sin3u>t  =  3/4  sin  ut  -  1/4  sin  Zu>t  (B.4) 

So  after  neglecting  higher  harmonics  one  gets 

XY  =  [  (XY) o  +  1/2  {XY),2  ]  +  [  (AT),  +  3/4  (XY),3 }  sin  u ,t  (B.5) 
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Appendix  C 

Multiplication  of  Three  Harmonic 
Quantities 


Some  of  the  governing  equations  contain  products  of  three  harmonic  quantities. 
Multiplication  of  three  harmonics  X,Y,Z  can  then  be  performed  as  a  series  of  two 
multiplications  involving  two  harmonic  quantities  as  follows. 

XYZ  =  ( XY)Z 

=  [(XY)0  +  (XY),  sin  u/t  +  (XY)j2  sin2  wt  +  (XY),3  sin3  ut] 

.  (Zo  +  Z ,  sinu>£  4-  Z, 2  sin2  uit  +  Z,3  sin3  ujt)  (C.l) 

=  (XYZ) o  +  (XYZ),  sin  u/i  +  (XYZ),2  sin2  ut  +  (XYZ),3  sin3 

where 


(XYZ)  o  =  (XY)oZo 
(XYZ).  =  (XY)0Z,  +  (XY),Zo 
(XYZU  =  (XY)0Z.2  +  (XY).Z,  +  (XY),2Z0 
(XYZ).  3  =  (XY)0Z.3  +  (XY).Z,2  +  (XY),2Zs  +  (XY),3Zo 


and 


(XY)0  =  X0Yo 

(XY).  =  XoY.  +  X.Yo 


37 


(XY),2  =  X0Y.2  +  X.Y.  +  Xa2Y0 

(xy).  3  =  x0y.3  +  jr.y.2  +  xi2ys  +  x.3y0 

as  before.  Once  again,  neglecting  higher  harmonics  one  gets 

XY  Z  =  [(XYZ)Q  +  l/2{XYZ),2}  +  [(XYZ),  +  3/4(XYZ),3}smut  (C.2) 
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!  Table  1:  Beam  Material  Properties  (AS4/3501-6) 

0/90)3,  Laminate 

t  =  1.49  x  10~3  m 

m  =  0.0683  kg/m 

Ip  =  5.13  x  10-6  kg.  m 

En  =  3.68  x  106  N 

Eu  =  0.26  x  106  N 

£33  =  2.9  x  105  N 

£44  =  0.183  N.  m2 

Ess  =  0.707  N.  m2 

£66  =  276  N.  m2 

[45/01,  Laminate 

t  =  1.49  x  10"3  m 

i  m  =  0.0238  kg/m 

Ip  =  1.66  x  10~6  kg.  m 

En  =  1-32  x  106  N 

£22  =  0.27  x  IQ6  N 

£33  =  1.0  x  105  N 

£44  =  0.0195  N.  m2 

Ess  =  0.0143  N.  m2 

£66  =  99.1  N.  m2 

En  =  1.0  x  106  N 

£4S  =  0.00632  N.  m2 

Note:  in  more  conventional  terms, 


E\\  —  EA  E22  —  GAn  £3 3  —  GAq 
E44  ~  GJ  E$s  —  ££)  -£'66  —  £7c 


£12  —  Extension-shear  coupling 

£14  ~  Extension-twist  coupling 

£4 5  ~  Bending-twist  coupling 
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Figure  5:  Second  Bending  Mode;[0/90)3,,0  mm  tip  deflection, Zs=10  mm 
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Figure  6:  Second  Bending  Mode;[0/90]3,,0  mm  tip  deflection, Zs= 100  mm 
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Figure  7:  First  Fore-and-Aft  Mode;[0/90|3,,0  mm  tip  deflection, Ys  =  •  rr.m 
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Figure  8:  First  Fore-and-Aft  Mode;,0  90_ 3, ,0  mm  tip  deflection, Ys=38  mm 
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Figure  11:  First  Bending  Mode;i0/90’3>,59  mm  tip  deflection, Zs=10 
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Figure  12:  First  Bending  Mode;[0/90]3,,59  mm  tip  deflection, Zs=200  mm 
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Figure  19:  First  Bending  Mode; [0 / 90]3» ,210  mm  tip  deflection, Zs=10  mm 
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Figure  21:  Second  Bending  Mode;[0/90]3,,210  mm  tip  deflection, Zs=  10  mm 
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Figure  22:  Second  Bending  Uode;[0/90]3,,210  mm  tip  deflection, Z*=  48  mm 
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Figure  27:  Frequency  vs.  Amplitude;  0/90)3,, 0  mm  tip  deflection 
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Figure  29:  Frequency  vs.  Amplitude;  0  90  3.. 210  mm  tip  deflection 
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Figure  32:  Frequency  vs.  Amplitude  w/o  and  w/  2nd  harmonics;  [O/90I3 
and  210  mm  tip  deflection 
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Figure  33:  Natural  Frequencies  of  r0/90]3.  Beam  as  a  Function  of  Tip  Deflection  (from 
Ref.  1) 
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Figure  34:  First  Bending  Mode;[45/0),,0  mm  tip  deflection, Zs=  10 
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Figure  35:  First  Bending  Mode;[45/0],,0  mm  tip  deflection, Zs=200  mm 


Figure  38:  First  Fore-and-Aft  Mode;[45/0j,,0  mm  tip  deflection, Ys=  1  n 
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Figure  39:  First  Fore-and-Aft  Mode;;45/0l,,0  mm  tip  deflection, Ys=  3.5 
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Figure  42:  First  Bending  Mode;j45/'0},,70  mm  tip  deflection, Zs=10  mm 
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Figure  43:  First  Bending  Mode;[45/0j,,70  mm  tip  deflection, Zs=200  mm 
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Figure  44:  Second  Bending  Mode;[45/0],,70  mm  tip  deflection, Zs= 10  mm 
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Figure  45:  Second  Bending  Mode;'45/0],,70  mm  tip  deflection, Zs=  70  mm 
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Figure  46:  First  Fore-and- Aft  Mode;[45/0],,70  mm  tip  deflection, Ys  =  10  mm 
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Figure  53:  Second  Bending  Mode;  45/0', ,203  mm  tip  deflection, Zs-- 45  mm 
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Figure  58:  Frequency  vs.  Aroplitude;;45/0l„0  mm  tip  deflection 
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Figure  62:  Centershift  v3.  Amplitude:;45/0!,,203  mm  tip  deflection 
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Figure  63:  Frequency  vs.  Amplitude  w  0  and  w  2nd  harmonics;  45  0  ,,70  mm 
and  203  mm  tip  deflection 
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